{
 "cells": [
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "# 遥感原理与应用\n",
    "\n",
    "## 1、DN\n",
    "\n",
    "**DN值（Digital Number，数字量化值）** 是遥感影像像元（像素）的原始数字记录值，表示传感器接收到的电磁波辐射能量的相对强度。\n",
    "\n",
    "### 1. DN值的本质\n",
    "\n",
    "- **传感器直接输出**：DN值是卫星/航空传感器记录的原始数字信号，未经物理量转换。\n",
    "- **整数存储**：通常以8-bit（0-255）、16-bit（0-65535）等整型格式存储。\n",
    "- **相对值**：仅代表辐射能量的相对强弱，无明确物理单位。\n",
    "\n",
    "### 2. DN值与物理量的关系\n",
    "\n",
    "需通过**辐射定标**转换为具有物理意义的量：\n",
    "\n",
    "| 转换目标                              | 公式                                                | 用途         |\n",
    "| :------------------------------------ | :-------------------------------------------------- | :----------- |\n",
    "| **辐射亮度** (Radiance)               | $L = DN \\times \\text{Gain} + \\text{Bias}$           | 大气校正输入 |\n",
    "| **反射率** (Reflectance)              |                                                     | 地表分析     |\n",
    "| **亮度温度** (Brightness Temperature) | $T = \\frac{\\ln\\left(\\frac{L}{K_1} + 1\\right)}{K_2}$ | 热红外波段   |\n",
    "\n",
    "> （$E_{\\text{sun}}$: 大气顶层太阳辐照度；θ: 太阳高度角；d: 日地距离修正系数）\n",
    "\n",
    "### 3. DN值的典型特征\n",
    "\n",
    "- **传感器依赖性**：不同传感器（如Landsat、Sentinel）的DN范围不同：\n",
    "  - Landsat 8/9 OLI: 0-65535（16-bit）\n",
    "  - Landsat 7 ETM+: 0-255（8-bit）\n",
    "- **波段差异性**：同一场景不同波段的DN值分布不同（如近红外DN通常高于短波红外）。\n",
    "\n",
    "## 2、Landsat\n",
    "\n",
    "**Landsat（陆地卫星）** 是由美国NASA和USGS联合运营的全球最长寿的地球观测卫星计划，自1972年发射第一颗卫星（Landsat 1）以来，持续提供中分辨率遥感影像，是遥感领域的基石。\n",
    "\n",
    "### 1. 核心特征\n",
    "\n",
    "| 特性         | 说明                                                         |\n",
    "| :----------- | :----------------------------------------------------------- |\n",
    "| **分辨率**   | 空间分辨率：15m（全色）-30m（多光谱） 时间分辨率：16天（重访周期） |\n",
    "| **光谱波段** | 8-11个波段（可见光-热红外）                                  |\n",
    "| **覆盖范围** | 全球陆地表面，幅宽185km×185km                                |\n",
    "| **数据级别** | L1（辐射定标）和L2（地表反射率/温度产品）                    |\n",
    "\n",
    "### 2. 主要卫星系列\n",
    "\n",
    "| 卫星            | 服役时间     | 关键升级                                 |\n",
    "| :-------------- | :----------- | :--------------------------------------- |\n",
    "| **Landsat 1-3** | 1972-1983    | 首开民用对地观测先河                     |\n",
    "| **Landsat 4-5** | 1982-2013    | 新增短波红外和热红外波段                 |\n",
    "| **Landsat 7**   | 1999-今      | ETM+传感器（全色15m）                    |\n",
    "| **Landsat 8-9** | 2013/2021-今 | OLI/TIRS传感器，新增海岸气溶胶和卷云波段 |\n",
    "\n",
    "### 3. 典型应用场景\n",
    "\n",
    "- **地表覆盖变化**：森林砍伐、城市扩张监测（如1990-2020年亚马逊雨林变化）\n",
    "- **农业**：作物长势评估（NDVI指数）、灌溉管理\n",
    "- **水资源**：湖泊面积变化（如咸海萎缩）\n",
    "- **灾害评估**：火灾烧伤严重度分析（dNBR指数）\n",
    "- **气候变化**：冰川退缩监测（如喜马拉雅冰川）\n",
    "\n",
    "### 4. 数据获取与处理\n",
    "\n",
    "- [USGS EarthExplorer](https://earthexplorer.usgs.gov/)\n",
    "- [Google Earth Engine](https://earthexplorer.usgs.gov/)\n",
    "- [NASA Earthdata](https://earthdata.nasa.gov/)\n",
    "- [地理空间数据云](https://www.gscloud.cn/home)\n",
    "\n",
    "### 5. 技术参数对比（Landsat 8/9 OLI）\n",
    "\n",
    "| 波段  | 波长范围(μm) | 主要用途           |\n",
    "| :---- | :----------- | :----------------- |\n",
    "| 1     | 0.43-0.45    | 海岸/气溶胶        |\n",
    "| 2     | 0.45-0.51    | 蓝（水体穿透）     |\n",
    "| 3     | 0.53-0.59    | 绿（植被活力）     |\n",
    "| 4     | 0.64-0.67    | 红（叶绿素吸收）   |\n",
    "| 5     | 0.85-0.88    | 近红外（生物量）   |\n",
    "| 6     | 1.57-1.65    | 短波红外1（水分）  |\n",
    "| 7     | 2.11-2.29    | 短波红外2（矿物）  |\n",
    "| 8     | 0.50-0.68    | 全色（15m）        |\n",
    "| 9     | 1.36-1.38    | 卷云检测           |\n",
    "| 10-11 | 10.6-12.5    | 热红外（地表温度） |\n",
    "\n",
    "### 6. 与其他卫星对比\n",
    "\n",
    "| 特性     | Landsat  | Sentinel-2 | MODIS     |\n",
    "| :------- | :------- | :--------- | :-------- |\n",
    "| 分辨率   | 30m      | 10-60m     | 250-1000m |\n",
    "| 重访周期 | 16天     | 5天        | 1-2天     |\n",
    "| 光谱波段 | 11个     | 13个       | 36个      |\n",
    "| 数据政策 | 完全开放 | 开放       | 开放      |\n",
    "\n"
   ],
   "id": "b419392d60192d2c"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "# Landsat数据预处理完整流程详解\n",
    "\n",
    "**数据可视化**:\n",
    "\n",
    "代码中包含了多种可视化方法，包括：\n",
    "\n",
    "- 单波段灰度图显示\n",
    "- RGB真彩色合成\n",
    "- 假彩色合成\n",
    "- 直方图分析\n",
    "\n",
    "**关键Python知识点**:\n",
    "\n",
    "1. **rasterio库**：用于处理栅格地理数据的主要库，可以读写GeoTIFF等格式\n",
    "2. **numpy数组操作**：处理多维数组数据（波段×高度×宽度）\n",
    "3. **geopandas**：处理矢量边界数据\n",
    "4. **matplotlib**：数据可视化\n",
    "5. **文件路径处理**：使用os.path处理文件路径\n",
    "6. **元数据解析**：从MTL文本文件中提取关键参数\n",
    "\n",
    "**预处理流程**:\n",
    "\n",
    "1. 读取原始单波段数据\n",
    "2. 将单波段堆叠成多波段\n",
    "3. 辐射校正（DN值→反射率）\n",
    "4. 使用边界裁剪研究区\n",
    "5. 保存预处理结果\n",
    "6. 可视化验证\n",
    "\n"
   ],
   "id": "8fdb12849bc5fac2"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## 1. 环境准备与数据理解\n",
    "\n",
    "- 数据文件分析:\n",
    "    - **1992年数据**：包含7个波段(B1-B7)的TIFF文件，以及关键的MTL.txt元数据文件\n",
    "    - **2021年数据**：包含多个波段的TIFF文件(应该是B1-B11)\n",
    "    - **新都区边界**：xinDu_2021_shp.shp等shapefile相关文件"
   ],
   "id": "4438bb5764c96733"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 1.1查看和可视化 Landsat 数据",
   "id": "9a92b71b5696d89f"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "import rasterio\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# 打开单个波段\n",
    "with rasterio.open('../data/landsat/LT51290391992229BJC00/LT51290391992229BJC00_B4.TIF') as src:\n",
    "    # 读取数据\n",
    "    data = src.read(1)  # 读取第一个波段\n",
    "\n",
    "    # 创建图形\n",
    "    plt.figure(figsize=(10, 10))\n",
    "\n",
    "    # 使用imshow显示并保存返回的图像对象\n",
    "    img = plt.imshow(data, cmap='gray')\n",
    "\n",
    "    plt.title('Band 4 (Near Infrared)')\n",
    "\n",
    "    # 添加colorbar，使用img作为mappable\n",
    "    plt.colorbar(img, label='反射率值')\n",
    "\n",
    "    plt.show()"
   ],
   "id": "aa7a53cff3b7ad18",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 解决中文乱码的问题\n",
    "plt.rcParams['font.sans-serif'] = ['SimHei']  # 设置黑体（Windows常用）\n",
    "plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题"
   ],
   "id": "d1ef608ddb67a8c4",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## 2. 1992年Landsat 4-5 TM数据预处理\n",
    "\n",
    "### 2.1 数据读取与波段堆叠\n",
    "1. **波段堆叠原理**：\n",
    "   - Landsat每个波段单独存储为一个TIFF文件\n",
    "   - 需要将多个单波段文件合并为多波段数组才能进行后续分析（如NDVI计算、分类等）\n",
    "2. **内存效率**：\n",
    "   - 使用`np.zeros`预分配内存避免重复扩容\n",
    "   - 指定`dtype=np.float32`保证数据精度同时节省内存\n",
    "3. **元数据保留**：\n",
    "   - 拷贝第一个波段的元数据并更新波段数，确保空间参考信息不丢失\n",
    "4. **扩展性**：\n",
    "   - `read_and_stack_bands`函数可复用处理其他年份或其他传感器的数据（只需修改输入路径）"
   ],
   "id": "e7b404993fc354a7"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "import rasterio\n",
    "from rasterio.mask import mask  # 明确导入mask函数\n",
    "import geopandas as gpd\n",
    "import os"
   ],
   "id": "4bf7c4b1513aa04a",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 设置1992年数据路径\n",
    "data_dir_1992 = \"../data/landsat/LT51290391992229BJC00\"  # Landsat文件路径\n",
    "# 使用列表推导式批量生成文件路径（避免手动写7行路径）\n",
    "bands_1992 = [os.path.join(data_dir_1992, f\"LT51290391992229BJC00_B{i}.TIF\") for i in range(1, 8)]\n",
    "\n",
    "\n",
    "def read_and_stack_bands(band_files):\n",
    "    \"\"\"\n",
    "    读取多个波段并堆叠成一个多波段图像\n",
    "    参数:\n",
    "        band_files: 波段文件路径列表\n",
    "    返回:\n",
    "        stacked: 堆叠后的多维数组(波段,高,宽)\n",
    "        meta: 元数据字典\n",
    "    \"\"\"\n",
    "    # 读取第一个波段获取元数据\n",
    "    with rasterio.open(band_files[0]) as src:\n",
    "        meta = src.meta\n",
    "        height, width = src.shape\n",
    "\n",
    "    # 初始化数组\n",
    "    stacked = np.zeros((len(band_files), height, width), dtype=np.float32)\n",
    "\n",
    "    # 读取所有波段\n",
    "    for i, band_file in enumerate(band_files):\n",
    "        with rasterio.open(band_file) as src:\n",
    "            stacked[i] = src.read(1)\n",
    "\n",
    "    # 更新元数据\n",
    "    meta.update(count=len(band_files))\n",
    "\n",
    "    return stacked, meta\n",
    "\n",
    "\n",
    "# 读取并堆叠1992年7个波段\n",
    "stacked_1992, meta_1992 = read_and_stack_bands(bands_1992)\n",
    "print(f\"1992年数据形状: {stacked_1992.shape} (波段, 高度, 宽度)\")"
   ],
   "id": "dbaa4e5f04165636",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "### 2.2 辐射校正（DN值转反射率）\n",
    "- 需要从MTL.txt文件中获取辐射校正参数。\n",
    "\n",
    "- Landsat数据的辐射校正，将原始的DN值（Digital Number，数字量化值）转换为地表反射率（Surface Reflectance）。\n",
    "- 辐射校正的公式为:\n",
    "    $\\text{地表反射率} = \\frac{\\pi \\times \\text{辐射亮度} \\times 100}{\\text{ESUN} \\times \\sin(\\text{太阳高度角})}$\n",
    "\n",
    "* 为什么需要辐射校正？\n",
    "  1. **DN值本身无物理意义**，只是传感器记录的原始数字\n",
    "  2. **反射率具有可比性**：\n",
    "     - 消除太阳高度角、大气条件的影响\n",
    "     - 使不同时间/地点的影像可以直接比较\n",
    "  3. **支持定量分析**：\n",
    "     - 植被指数计算（如NDVI）\n",
    "     - 地表覆盖分类"
   ],
   "id": "e413a500d81f010e"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def parse_mtl(mtl_file):\n",
    "    \"\"\"\n",
    "    解析Landsat MTL文件获取辐射校正参数\n",
    "    参数:\n",
    "        mtl_file: MTL文件路径\n",
    "    返回:\n",
    "        包含关键参数的字典\n",
    "    \"\"\"\n",
    "    params = {}\n",
    "    with open(mtl_file, 'r') as f:\n",
    "        for line in f:\n",
    "            if '=' in line:\n",
    "                key, value = line.split('=', 1)\n",
    "                key = key.strip()\n",
    "                value = value.strip().strip('\"')\n",
    "                params[key] = value\n",
    "    return params\n",
    "\n",
    "\n",
    "# 解析MTL文件\n",
    "mtl_1992 = parse_mtl(os.path.join(data_dir_1992, \"LT51290391992229BJC00_MTL.txt\"))\n",
    "\n",
    "# 获取太阳高度角(度)\n",
    "sun_elevation = float(mtl_1992['SUN_ELEVATION'])\n",
    "\n",
    "# Landsat5 TM各波段的太阳辐照度(ESUN) - 单位: W/(m²·μm)\n",
    "esun_tm = [1957.0, 1826.0, 1554.0, 1036.0, 215.0, 80.67, 1362.0]\n",
    "\n",
    "\n",
    "def dn_to_reflectance(band_array, band_index, sun_elevation, esun):\n",
    "    \"\"\"\n",
    "    将DN值转换为地表反射率\n",
    "    参数:\n",
    "        band_array: 波段数组\n",
    "        band_index: 波段索引(0-based)\n",
    "        sun_elevation: 太阳高度角(度)\n",
    "        esun: 太阳辐照度列表\n",
    "    返回:\n",
    "        反射率数组\n",
    "    \"\"\"\n",
    "    # 从MTL获取增益和偏置\n",
    "    gain = float(mtl_1992[f'RADIANCE_MULT_BAND_{band_index + 1}'])\n",
    "    bias = float(mtl_1992[f'RADIANCE_ADD_BAND_{band_index + 1}'])\n",
    "\n",
    "    # 计算辐射亮度\n",
    "    radiance = band_array * gain + bias\n",
    "\n",
    "    # 计算反射率\n",
    "    reflectance = (np.pi * radiance * 100) / (esun[band_index] * np.sin(np.radians(sun_elevation)))\n",
    "\n",
    "    return reflectance\n",
    "\n",
    "\n",
    "# 对每个波段进行辐射校正\n",
    "for i in range(stacked_1992.shape[0]):\n",
    "    stacked_1992[i] = dn_to_reflectance(stacked_1992[i], i, sun_elevation, esun_tm)"
   ],
   "id": "e8a267a9a241a226",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# mtl_1992\n",
    "stacked_1992"
   ],
   "id": "692a9763084cade5",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 2.3 使用shapefile裁剪研究区域\n",
   "id": "c987bfe02ae57842"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 读取新都区边界\n",
    "xindu_shp = gpd.read_file(\"../data/shape_file/xinDu_2021_shp/xinDu_2021_shp.shp\")\n",
    "\n",
    "\n",
    "def clip_raster_with_shapefile(raster, meta, shapefile):\n",
    "    \"\"\"\n",
    "    使用shapefile裁剪栅格数据\n",
    "    参数:\n",
    "        raster: 栅格数据数组 (形状: [波段数, 高度, 宽度])\n",
    "        meta: 元数据字典\n",
    "        shapefile: GeoDataFrame包含边界\n",
    "    返回:\n",
    "        裁剪后的栅格数据和更新后的元数据\n",
    "    \"\"\"\n",
    "    # 创建临时目录（如果不存在）\n",
    "    os.makedirs(\"../data/temp\", exist_ok=True)\n",
    "    temp_file = \"../data/temp/temp_stack.tif\"\n",
    "\n",
    "    try:\n",
    "        # 临时保存堆叠的栅格数据\n",
    "        with rasterio.open(temp_file, 'w', **meta) as dst:\n",
    "            dst.write(raster)\n",
    "\n",
    "        # 使用rasterio进行裁剪\n",
    "        with rasterio.open(temp_file) as src:\n",
    "            # 确保shapefile与栅格同坐标系\n",
    "            shapefile = shapefile.to_crs(src.crs)\n",
    "\n",
    "            # 获取几何图形（必须是闭合多边形）\n",
    "            geoms = shapefile.geometry.values  # 直接获取GeoSeries的值\n",
    "\n",
    "            # 执行裁剪（注意：mask函数需要几何图形列表）\n",
    "            out_image, out_transform = mask(src, geoms, crop=True, all_touched=True)\n",
    "            out_meta = src.meta.copy()\n",
    "\n",
    "        # 更新元数据\n",
    "        out_meta.update({\n",
    "            \"height\": out_image.shape[1],\n",
    "            \"width\": out_image.shape[2],\n",
    "            \"transform\": out_transform\n",
    "        })\n",
    "\n",
    "        return out_image, out_meta\n",
    "\n",
    "    finally:\n",
    "        # 确保临时文件被删除\n",
    "        if os.path.exists(temp_file):\n",
    "            os.remove(temp_file)\n",
    "\n",
    "\n",
    "# 裁剪1992年数据\n",
    "clipped_1992, meta_1992 = clip_raster_with_shapefile(stacked_1992, meta_1992, xindu_shp)"
   ],
   "id": "858760d5b56cfd7a",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "#### 查看裁剪后的栅格数据结果\n",
    "\n",
    "1. 使用matplotlib快速查看单个波段\n",
    "2. 查看RGB真彩色合成（波段3/2/1组合）\n",
    "3. 检查基本统计信息"
   ],
   "id": "dea470f6528fc89a"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# 显示第一个波段（索引0）\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.imshow(clipped_1992[0], cmap='gray')  # 可以修改索引查看不同波段\n",
    "plt.title('裁剪后的波段1 (1992年)')\n",
    "plt.colorbar(label='反射率值')\n",
    "plt.show()"
   ],
   "id": "e1dcdfb024fc7056",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 选择波段 (注意波段顺序是BGR)\n",
    "rgb_bands = [2, 1, 0]  # Landsat TM的波段3(R)、2(G)、1(B)\n",
    "\n",
    "# 提取并标准化波段\n",
    "rgb = np.stack([\n",
    "    clipped_1992[rgb_bands[0]],  # 红\n",
    "    clipped_1992[rgb_bands[1]],  # 绿\n",
    "    clipped_1992[rgb_bands[2]]  # 蓝\n",
    "], axis=-1)\n",
    "\n",
    "# 标准化到0-1范围\n",
    "rgb_norm = (rgb - rgb.min()) / (rgb.max() - rgb.min())\n",
    "\n",
    "# 显示\n",
    "plt.figure(figsize=(12, 10))\n",
    "plt.imshow(rgb_norm)\n",
    "plt.title('1992年真彩色合成 (RGB)')\n",
    "plt.axis('off')\n",
    "plt.show()"
   ],
   "id": "7f3942aa7f7ab707",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print(\"裁剪结果信息:\")\n",
    "print(f\"波段数: {clipped_1992.shape[0]}\")\n",
    "print(f\"空间范围: {meta_1992['transform'] * (0, 0)} 到 \"\n",
    "      f\"{meta_1992['transform'] * (meta_1992['width'], meta_1992['height'])}\")\n",
    "print(\"各波段值范围:\")\n",
    "for i in range(clipped_1992.shape[0]):\n",
    "    print(f\"波段{i + 1}: {np.nanmin(clipped_1992[i]):.2f} ~ {np.nanmax(clipped_1992[i]):.2f}\")"
   ],
   "id": "7f93872dc22a7c51",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 2.4 保存预处理结果",
   "id": "e68687dcf387aec0"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 更新元数据以反映反射率值范围\n",
    "meta_1992.update(dtype='float32', nodata=0)\n",
    "\n",
    "# 保存预处理后的1992年数据\n",
    "output_1992 = \"../data/output_preprocessed/preprocessed_landsat5_1992.tif\"\n",
    "with rasterio.open(output_1992, 'w', **meta_1992) as dst:\n",
    "    dst.write(clipped_1992)\n",
    "\n",
    "print(f\"1992年数据预处理完成，已保存到: {output_1992}\")"
   ],
   "id": "75541b6ebcf436de",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "#### 可视化预处理后的TIFF文件\n",
    "1. 单波段灰度图\n",
    "2. 真彩色合成（RGB）"
   ],
   "id": "2272f32e1ed2e1ac"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "import rasterio\n",
    "from rasterio.plot import show\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "with rasterio.open(output_1992) as src:\n",
    "    # 读取第1个波段（注意：波段索引从1开始）\n",
    "    band_data = src.read(1)  # 读取第一个波段\n",
    "\n",
    "    # 创建图形\n",
    "    plt.figure(figsize=(10, 8))\n",
    "\n",
    "    # 使用imshow显示数据\n",
    "    img = plt.imshow(band_data, cmap='gray')\n",
    "\n",
    "    # 添加颜色条\n",
    "    plt.colorbar(img, label='反射率值 (×100)')\n",
    "    plt.title('1992年预处理数据 - 波段1')\n",
    "    plt.show()"
   ],
   "id": "aefd14d9ef509686",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "with rasterio.open(output_1992) as src:\n",
    "    # Landsat 5 TM的波段顺序：B1(蓝), B2(绿), B3(红)\n",
    "    rgb = src.read([3, 2, 1])  # 读取为[R,G,B]顺序\n",
    "\n",
    "    # 标准化到0-1范围\n",
    "    rgb_norm = (rgb - rgb.min()) / (rgb.max() - rgb.min())\n",
    "\n",
    "    plt.figure(figsize=(12, 10))\n",
    "    plt.imshow(rgb_norm.transpose(1, 2, 0))  # 转为[高,宽,波段]顺序\n",
    "    plt.title('1992年真彩色合成 (波段3/2/1)')\n",
    "    plt.axis('off')\n",
    "    plt.show()"
   ],
   "id": "8e2691d4de441a47",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## 3. 2021年Landsat 8数据预处理",
   "id": "a55dd11147797957"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 3.1 路径设置与验证",
   "id": "d1d04c26e69ef90b"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "data_dir_2021 = \"../data/landsat/LC81290392021036LGN00\"\n",
    "bands_2021 = [os.path.join(data_dir_2021, f\"LC08_L1TP_129039_20210205_20210304_01_T1_B{i}.TIF\")\n",
    "              for i in [1, 2, 3, 4, 5, 6, 7, 9]]  # Landsat8 OLI波段\n",
    "\n",
    "# 验证文件存在性\n",
    "print(f\"第一个波段文件存在: {os.path.exists(bands_2021[0])}\")\n"
   ],
   "id": "289976be77c00652",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 3.2 MTL解析",
   "id": "ca81c454671ba82"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def parse_mtl_2021(mtl_file):\n",
    "    \"\"\"增强版MTL解析，处理科学计数法和异常值\"\"\"\n",
    "    params = {}\n",
    "    with open(mtl_file, 'r') as f:\n",
    "        for line in f:\n",
    "            if '=' in line:\n",
    "                key, value = line.split('=', 1)\n",
    "                key = key.strip()\n",
    "                value = value.strip().strip('\"')\n",
    "                try:\n",
    "                    # 特殊处理科学计数法\n",
    "                    if 'e' in value.lower():\n",
    "                        params[key] = float(value)\n",
    "                    else:\n",
    "                        # 尝试转换为int或float\n",
    "                        params[key] = float(value) if '.' in value else int(value)\n",
    "                except (ValueError, TypeError):\n",
    "                    params[key] = value\n",
    "    return params\n"
   ],
   "id": "908455a134bcd698",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 3.3 数据读取与堆叠",
   "id": "aa31fa4ec0ba34bc"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def read_and_stack_bands(band_files):\n",
    "    \"\"\"读取并堆叠多个波段\"\"\"\n",
    "    with rasterio.open(band_files[0]) as src:\n",
    "        meta = src.meta\n",
    "        height, width = src.shape\n",
    "\n",
    "    # 初始化数组时使用float32防止溢出\n",
    "    stacked = np.zeros((len(band_files), height, width), dtype=np.float32)\n",
    "\n",
    "    for i, band_file in enumerate(band_files):\n",
    "        with rasterio.open(band_file) as src:\n",
    "            stacked[i] = src.read(1).astype(np.float32)  # 显式转换为float32\n",
    "\n",
    "    meta.update(count=len(band_files), dtype='float32')\n",
    "    return stacked, meta\n",
    "\n",
    "\n",
    "# 执行读取\n",
    "stacked_2021, meta_2021 = read_and_stack_bands(bands_2021)\n",
    "print(f\"\\n原始DN值范围: {stacked_2021.min()}~{stacked_2021.max()}\")\n"
   ],
   "id": "9c82112e69b6ba45",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 3.4 辐射校正",
   "id": "213c343a7ad462a7"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "mtl_path = os.path.join(data_dir_2021, \"LC08_L1TP_129039_20210205_20210304_01_T1_MTL.txt\")\n",
    "mtl_2021 = parse_mtl_2021(mtl_path)\n",
    "\n",
    "# 验证关键参数\n",
    "print(\"\\nMTL参数验证:\")\n",
    "print(\"RADIANCE_MULT样例:\", {k: v for k, v in mtl_2021.items() if \"RADIANCE_MULT\" in k})\n",
    "print(\"RADIANCE_ADD样例:\", {k: v for k, v in mtl_2021.items() if \"RADIANCE_ADD\" in k})\n",
    "\n",
    "\n",
    "def dn_to_radiance(band_array, band_index):\n",
    "    \"\"\"DN转辐射亮度（更安全的实现）\"\"\"\n",
    "    band_num = [1, 2, 3, 4, 5, 6, 7, 9][band_index]\n",
    "    try:\n",
    "        mult = mtl_2021[f'RADIANCE_MULT_BAND_{band_num}']\n",
    "        add = mtl_2021[f'RADIANCE_ADD_BAND_{band_num}']\n",
    "        radiance = band_array * mult + add\n",
    "        return np.clip(radiance, 0, None)  # 剔除负值\n",
    "    except KeyError:\n",
    "        print(f\"警告: 波段{band_num}参数缺失，使用默认转换\")\n",
    "        return band_array * 0.1  # 应急默认值\n",
    "\n",
    "\n",
    "# 执行辐射校正\n",
    "print(\"\\n辐射校正过程:\")\n",
    "for i in range(stacked_2021.shape[0]):\n",
    "    band_idx = [1, 2, 3, 4, 5, 6, 7, 9][i]\n",
    "    print(f\"波段{band_idx}校正前: {stacked_2021[i].min():.2f}~{stacked_2021[i].max():.2f}\")\n",
    "    stacked_2021[i] = dn_to_radiance(stacked_2021[i], i)\n",
    "    print(f\"波段{band_idx}校正后: {stacked_2021[i].min():.4f}~{stacked_2021[i].max():.4f}\")\n"
   ],
   "id": "b97e69e283def981",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 3.5 数据裁剪",
   "id": "b37252720548a1a6"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 裁剪2021年数据\n",
    "clipped_2021, meta_2021 = clip_raster_with_shapefile(stacked_2021, meta_2021, xindu_shp)\n"
   ],
   "id": "b6606a5d6bbafb2",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 3.6保存结果",
   "id": "8db4a3e5609f0c8c"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "output_dir = \"../data/output_preprocessed\"\n",
    "os.makedirs(output_dir, exist_ok=True)\n",
    "# output_2021 = os.path.join(output_dir, \"landsat8_2021_radiance.tif\")\n",
    "output_2021 = os.path.join(output_dir, \"preprocessed_landsat8_2021.tif\")\n",
    "\n",
    "# 更新元数据\n",
    "meta_2021.update(dtype='float32', nodata=0)\n",
    "\n",
    "with rasterio.open(output_2021, 'w', **meta_2021) as dst:\n",
    "    dst.write(clipped_2021)\n",
    "print(f\"\\n辐射亮度数据已保存至: {output_2021}\")"
   ],
   "id": "c72df6da602daccd",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 3.7 结果验证",
   "id": "c925264b7665298d"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def plot_band_histogram(data, band_name):\n",
    "    \"\"\"绘制波段直方图\"\"\"\n",
    "    plt.figure()\n",
    "    plt.hist(data.flatten(), bins=50, range=(0, np.percentile(data, 99)))\n",
    "    plt.title(f\"{band_name}值分布\")\n",
    "    plt.xlabel(\"辐射亮度值\")\n",
    "    plt.ylabel(\"像素数量\")\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "print(\"\\n结果验证:\")\n",
    "with rasterio.open(output_2021) as src:\n",
    "    # 打印统计信息\n",
    "    for i in range(1, src.count + 1):\n",
    "        band = src.read(i)\n",
    "        print(f\"波段{i}: {band.min():.4f}~{band.max():.4f} (均值:{band.mean():.4f})\")\n",
    "\n",
    "    # 可视化近红外波段（B5）\n",
    "    band5 = src.read(5)\n",
    "    fig, ax = plt.subplots(figsize=(10, 8))\n",
    "\n",
    "    # 使用imshow替代show()，并保存返回的图像对象\n",
    "    img = ax.imshow(band5, cmap='gray', vmax=np.percentile(band5, 98))\n",
    "    plt.title('近红外波段（B5）辐射亮度')\n",
    "\n",
    "    # 使用fig.colorbar()并传入img对象\n",
    "    fig.colorbar(img, ax=ax, label='辐射亮度 (W/(m²·sr·μm))')\n",
    "    plt.show()\n",
    "\n",
    "    # 绘制直方图\n",
    "    plot_band_histogram(band5, \"B5\")"
   ],
   "id": "3d599a7bea896e44",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "#### 可视化预处理后的TIFF文件\n",
    "1. 单波段灰度图\n",
    "2. 真彩色合成（RGB）"
   ],
   "id": "c501c8f3d6479e0e"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 读取 GeoTIFF 文件\n",
    "import rasterio\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# 读取保存的 GeoTIFF 文件\n",
    "output_2021 = \"../data/output_preprocessed/preprocessed_landsat8_2021.tif\"\n",
    "\n",
    "with rasterio.open(output_2021) as src:\n",
    "    data = src.read()  # 读取所有波段数据\n",
    "    profile = src.profile  # 获取元数据\n",
    "    bounds = src.bounds  # 获取地理范围\n",
    "\n",
    "print(\"数据形状 (波段数, 高度, 宽度):\", data.shape)\n",
    "print(\"元数据:\", profile)\n",
    "print(\"地理范围:\", bounds)"
   ],
   "id": "44ad3a12594fef42",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 可视化单波段（灰度图）\n",
    "# 假设第一个波段是可见光波段（如蓝光）\n",
    "plt.figure(figsize=(10, 8))\n",
    "plt.imshow(data[0], cmap='gray')  # 显示第一个波段\n",
    "plt.colorbar(label='辐射亮度值')\n",
    "plt.title(\"Landsat 8 2021 - 波段 1 (灰度图)\")\n",
    "plt.axis('off')\n",
    "plt.show()"
   ],
   "id": "ea82924169ee0465",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 假设波段顺序是 [蓝, 绿, 红, 近红外...]，我们取第 2、3、4 波段作为 RGB\n",
    "rgb_bands = data[[2, 1, 0]]  # 调整顺序为 R, G, B\n",
    "\n",
    "# 归一化到 0-1 范围以便显示\n",
    "rgb_normalized = (rgb_bands - np.min(rgb_bands)) / (np.max(rgb_bands) - np.min(rgb_bands))\n",
    "\n",
    "plt.figure(figsize=(10, 8))\n",
    "plt.imshow(np.transpose(rgb_normalized, (1, 2, 0)))  # 转换为 (高度, 宽度, 波段) 格式\n",
    "plt.title(\"Landsat 8 2021 - 真彩色 RGB\")\n",
    "plt.axis('off')\n",
    "plt.show()"
   ],
   "id": "51f882285045d097",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## 4、结果验证与可视化",
   "id": "6cdb7b2bf7a5a335"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 4.1 检查输出文件\n",
   "id": "6aa1b3c279237b47"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def inspect_output(file_path):\n",
    "    \"\"\"检查输出文件的基本信息\"\"\"\n",
    "    with rasterio.open(file_path) as src:\n",
    "        print(f\"\\n文件: {file_path}\")\n",
    "        print(f\"波段数: {src.count}\")\n",
    "        print(f\"空间分辨率: {src.transform[0]}米\")\n",
    "        print(f\"范围: {src.bounds}\")\n",
    "        print(f\"坐标系: {src.crs}\")\n",
    "\n",
    "\n",
    "inspect_output(output_1992)\n",
    "inspect_output(output_2021)"
   ],
   "id": "90f09b9599da1344",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 4.2 可视化波段组合",
   "id": "4032fcb4bf99150e"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# 可视化1992年假彩色合成(543波段)\n",
    "def plot_rgb(raster, bands, title):\n",
    "    \"\"\"显示RGB合成图像\"\"\"\n",
    "    rgb = np.dstack((raster[bands[0]], raster[bands[1]], raster[bands[2]]))\n",
    "    # 拉伸对比度\n",
    "    rgb = np.clip(rgb * 3, 0, 1)  # 调整系数以获得最佳显示效果\n",
    "\n",
    "    plt.figure(figsize=(10, 10))\n",
    "    plt.imshow(rgb)\n",
    "    plt.title(title)\n",
    "    plt.axis('off')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "# 1992年543波段(假彩色)\n",
    "with rasterio.open(output_1992) as src:\n",
    "    raster_1992 = src.read()\n",
    "plot_rgb(raster_1992, [4, 3, 2], \"1992年假彩色合成(543波段)\")\n",
    "\n",
    "# 2021年543波段(假彩色)\n",
    "with rasterio.open(output_2021) as src:\n",
    "    raster_2021 = src.read()\n",
    "plot_rgb(raster_2021, [4, 3, 2], \"2021年假彩色合成(543波段)\")"
   ],
   "id": "39392f06edde31bb",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "id": "5f3110392c7807a9",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "id": "9b39e9aac3e3a0a3",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "id": "688022e0c710881b",
   "outputs": [],
   "execution_count": null
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
